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Section  1  -  TFE  1-D  Thermal  Model 


1.1.  Introduction 


The  TOPAZ-II  space  nuclear  reactor,  Ya-21u,  was  tested  as  part  of  the  Topaz  International 
Program  in  Albuquerque,  NM  during  the  period  August  1993  to  March  1995.  It  was  also  tested 
in  Russia  prior  to  its  arrival  in  the  U.S.  Thermal  testing  was  carried  out  using  tungsten  electrical 
heaters  (TISA  heaters)  to  simulate  the  nuclear  fuel.  No  nuclear  testing  has  been  carried  out  on 
this  system.  Details  of  the  TOPAZ-II  reactor  are  given  by  the  U.S.  TOPAZ-II  Flight  Safety  Team 
(1992),  Wold  et.  al.  (1994)  and  Paramonov  and  El-Genk  (1994).  During  the  testing  of  Ya-21u, 
the  reactor  was  subjected  to  a  variety  of  extreme  conditions,  including  a  shake  and  vibration  test, 
air  leaks  into  the  TFEs  and  thermal  shocks.  Details  of  U.S.  testing  of  Ya-21u  are  given  by 
Schmidt  et.  al.  (1995).  Over  the  duration  of  these  tests  the  electrical  output  parameters  of  the 
reactor  changed  significantly,  presumably  as  a  result  of  changes  to  internal  TFE  parameters  such 
as  effective  emissivity  of  the  inter-electrode  gap,  emitter  and  collector  work  functions  and  gap 
dimensions.  The  individual  TFE  characteristics  depend  sensitively  on  these  parameters  and  the 
final  electrical  output  is  the  result  of  a  delicate  balance  between  a  number  of  competing  effects. 
In  order  to  better  understand  the  factors  resulting  in  the  output  degradation  it  proved  necessary  to 
develop  computer  models  of  the  various  internal  TFE  processes.  An  additional  need  arose  to 
determine  the  effect  of  changes  in  NaK  coolant  temperature,  such  as  occurred  when  the  reactor 
was  run  without  the  radiator  thermal  shield.  The  computer  codes  described  in  this  report  were 
developed  to  meet  these  needs.  Their  primary  function  is  to  model  heat  transfer  and  thermionic 
emission  within  a  single  TFE.  The  theoretical  treatment  is  exclusively  one  dimensional.  This 
limitation  is  not  thought  to  be  particularly  important  for  our  primary  applications.  There  is  no 
doubt  that  a  one  dimensional  treatment  is  capable  of  indicating  qualitative  trends  in  the  variation 
of  parameters,  and  the  quantitative  limitations  of  a  one  dimensional  treatment  (due  to  TFE  'end 
effects')  are  not  too  severe,  in  view  of  the  actual  TFE  geometry.  In  our  models  heat  transfer 
within  the  TFE  is  calculated  out  to  the  NaK  coolant  channel  and  the  temperature  here  is  used  as  a 
boundary  condition.  A  more  comprehensive  model  of  the  entire  reactor  system,  including  coolant 
flow  and  reactor  neutronics,  is  described  by  Paramonov  and  El-Genk  (1994).  However  this 
degree  of  sophistication  is  not  necessary  for  the  applications  described  above.  By  eliminating 
details  of  the  rest  of  the  reactor  systems  it  has  been  possible  to  work  with  more  detailed  and 
realistic  codes  for  generating  the  current-voltage  characteristics  of  TFEs  under  a  given  set  of 
operating  conditions.  This  in  turn  should  lead  to  more  reliable  estimates  of  internal  TFE 
parameters. 


1.2.  Theory 

The  heat  transfer  models  combine  the  thermionic  emission  process  with  the  heat  transfer 
properties  of  the  TFE.  They  calculate  the  output  current  density  (Jout)  given  the  input  thermal 
power  (qFuEL)>  output  volts  (Vout),  cesium  pressure  (PCs)  and  NaK  flow  rate  (m).  Unless  the 
reactor  is  running  at  very  low  current  densities,  the  heat  transfer  across  the  electrode  gap  as  a 
result  of  the  transfer  of  electrons  from  the  emitter  to  the  collector  is  a  significant  proportion  of 


the  total  heat  transfer..  The  electron  current  is  a  very  sensitive  function  of  temperature  in  the  TFE, 
especially  the  emitter  temperature.  Hence  an  iterative  procedure  is  required  to  calculate  the 
temperatures  of  the  components  of  the  TFE  and  the  output  current  for  a  given  set  of  operating 
conditions.  This  procedure  alternately  balances  the  heat  fluxes  into  and  out  of  each  electrode. 

Three  codes  are  described  in  the  following  work.  These  are  IV ID,  EQ1D  and  NAK1D.  IV ID  is  a 
program  for  calculating  I-V  characteristics  only,  given  the  emitter  and  collector  temperatures  (TE, 
Tc),  Pcs,  and  Vout.  It  has  the  option  to  allow  for  voltage  drops  in  the  electrodes,  due  to  resistive 
losses.  For  this  it  also  needs  a  value  for  the  NaK  temperature  (TNaK),  which  is  used  in  calculating 
the  resistivity  of  the  electrodes  along  their  entire  length.  It  does  not  consider  heat  transfer  nor 
calculate  equilibrium  temperatures.  EQ1D  is  used  in  calculating  equilibrium  values  of  TE,  Tc  and 
Jout  in  operating  conditions  similar  to  those  used  on  Ya-21u.  It  is  used  to  produce  a  calibration 
between  the  heat  flux  to  the  NaK  (QNaK)  and  TNaK,  from  experimental  data  of  the  relationship 
between  qFuEL  (f°r  1  TFE)  and  TNa[C  (see  section  1.2.2).  NAK1D  is  the  most  general  program, 
calculating  equilibrium  values  of  TE,  Tc,  TNaK  and  Jout  in  a  TOPAZ-II  TFE  under  any  physically 
reasonable  operating  conditions.  It  then  allows  the  user  to  change  TNaK  and  calculate  the  new 
equilibrium  conditions.  This  is  useful  for  considering  the  changes  which  occurred  when  Ya-21u 
was  run  without  the  thermal  shield,  resulting  in  NaK  temperatures  up  to  ~  60  K  lower  than  when 
the  thermal  shield  was  in  place.  Both  heat  transfer  models  use  the  same  IV  characteristic  code  as 
IV1D.  The  main  components  of  the  model  are  described  below  in  section  1.2.2. 

1.2.1  Units 

Unless  otherwise  stated,  the  following  units  have  been  used  in  this  work,  in  both  the  codes  and 
the  theory  described  below: 


Quantity 

Symbol 

Units 

Temperature 

T 

K 

Current  density 

J 

A  cm_2 

Current 

I 

A 

Voltage 

V 

V 

Work  function 

<i> 

eV 

Power 

q 

W 

Heat  flux 

Q 

W  cm'z 

Area 

A 

cm2 

Length 

l 

cm 

Radius,  diameter 

r,  d 

cm 

Gap  spacing 

d 

cm 

Pressure 

P 

Torn 

Thermal  conductivity 

k,  X 

Electrical  resistivity 

P 

Q  cm 

Specific  heat  capacity 

cP 

IQSjgQUB 

Mass  flow  rate 

m 

kg  s'1 

1.2.2  Heat  transfer 


The  heat  transfer  models  NAK1D  and  EQ1D  both  find  equilibrium  conditions  by  using  a 
Newton-Raphson  routine  to  adjust  the  temperature  of  each  electrode  until  the  heat  flux  into  the 
electrode  equals  that  out  of  the  electrode.  This  is  done  alternately  for  the  emitter  and  collector 
until  the  temperatures  and  output  current  converge.  The  heat  flux  (i.e.  power  per  unit  area)  from 
the  TIS  A  to  the  inner  surface  of  the  emitter  (QFUEL)  is  calculated  from  qFUEL.  The  active  length  of 
the  TFE,  over  which  heat  transfer  from  the  TISA  has  been  assumed  to  take  place,  /,  is  30  cm.  The 
loss  from  qFUEL  as  a  result  of  resistive  losses  in  the  leads  is  taken  as  12%  (Steppenov  1991).  The 
radii  of  the  TFE  components  are  shown  in  figure  1.1  (Paramonov  and  El-Genk  1994). 
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.980  cm 
.030  cm 
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.  1 90  cm 
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.295  cm 


(1.1) 

Three  mechanisms  of  heat  transfer  from  the  emitter  to  the  collector  are  significant.  From  largest 
to  smallest  these  are  radiative,  electron  transfer  and  conductive  heat  transfer  through  the  Cs 
vapour.  Since  the  heat  transfer  via  electron  transfer  is  significant  under  typical  TFE  operating 
conditions,  this  effectively  couples  the  heat  transfer  and  I-V  characteristic  parts  of  the  model. 
The  radiative  heat  transfer  Qrad  is  given  by 

.  0-2) 

where  a  is  the  Stefan-Boltzmann  constant  and  seff  is  the  effective  emissivity  of  the  inter-electrode 
gap.  This  depends  on  the  emissivities  of  both  the  collector  and  the  emitter.  The  values  used  for 
seff  are  discussed  in  section  1 .2.2. 1 . 
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Figure  1.1:  Radii  of  TOPAZ-II  TFE  components. 
Hence  QFUEL  is  given  by 
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The  heat  flux  from  the  emitter  due  to  electron  cooling,  QeE,  depends  on  the  mode  in  which  the 
TFE  is  operating  (McVey  and  Rasor,  Dahlberg  et.  al.  1994,  Rasor  Associates  1991).  QeE  is 
calculated  in  the  IV  characteristic  part  of  the  program  and  will  not  be  discussed  further  here. 

The  corresponding  heat  flux  due  to  electron  heating  of  the  collector,  QeC,  is  given  by 

QeC=QeE-Jo«Ku,  •  0-3) 

The  heat  flux  due  to  conduction  through  the  Cs  vapour,  QCs,  is  given  by  Hatsopoulos  and 
Gyftopoulos  (1973)  as 

n  _  ^w(^n>X^~^c)  / 1 

d  + 1.15  x  10"5[(r£  +  Tc)/pa]  •  1  '  ’ 

where  d  is  the  emitter-collector  spacing  in  cm,  and  kmw(Tmw)  is  the  thermal  conductivity  of  Cs. 
This  is  a  function  of  the  weighted  mean  temperature  Tmw,  given  by  Hatsopoulos  and  Gyftopoulos 
(1973)  as 


T, . = 


n-ji 


A 


TE-Z 


c  J 


(1.5) 


The  variation  of  Xmw  with  Tmvv  is  shown  in  figure  CIO,  together  with  the  analytical  expression 
representing  this  data. 

Heat  transfer  from  the  collector  to  the  NaK  is  dominated  by  conduction,  since  the  temperatures 
are  much  lower  than  at  the  emitter.  Only  conduction  has  been  considered  in  these  models. 
Conduction  through  the  plasma-sprayed  A1203  collector  coating,  the  He  gap  and  the  stainless 
steel  inner  wall  of  the  NaK  channel  have  been  considered.  Finally,  the  convective  heat  transfer  to 
the  NaK  must  be  included.  Of  these  contributions,  the  He  gap  provides  the  greatest  thermal 
resistance.  Conduction  through  this  layer  structure  is  best  described  by  combining  the  thermal 
conductivities  of  each  layer  into  one  parameter.  Hence  an  effective  thermal  conductivity  per  unit 
length  of  the  active  zone  of  the  TFE,  Kmean,  is  defined  such  that 


Qncik  ~  Zmean  (Tc ,  TNaK  )(TC  TJvqk) 


(1.6) 


Kmean  is  a  function  of  both  Tc  and  TnoK’  since  the  thermal  conductivities  of  the  components 
between  the  collector  and  NaK  channel  are  all  functions  of  temperature.  Each  of  the  components 
can  be  thought  of  as  a  thermal  resistance.  For  conduction  in  the  alumina,  He  and  stainless  steel, 
the  resistance  of  that  element,  R^,  is  given  by  Kreith  and  Black  (1980)  as 

n  ln(r„+1/r„) 


2ukJ 


(1.7) 


where  /  is  the  length  over  which  conduction  takes  place  (i.e.  the  active  zone  of  the  TFE)  and  k,,  is 
the  thermal  conductivity  of  that  material.  With  reference  to  with  figure  1.1,  the  value  of  n  is  4  for 
the  alumina,  5  for  the  He  and  6  for  the  stainless  steel  layer.  Plots  of  the  thermal  conductivity  for 
alumina  (Paramonov  and  El-Genk  1994),  He  (CRC  1976)  and  stainless  steel  (Grigorieva  and 
Melikhova  1991),  as  a  function  of  temperature,  are  shown  in  figures  C5,  C6  and  C7  respectively. 
The  analytical  expressions  representing  these  data,  which  have  been  used  in  the  codes,  are  also 
shown.  The  thermal  resistance  of  convective  heat  transfer  to  the  NaK,  R7,  is  given  by  Kreith  and 
Black  (1980)  as 


*7  = 


1 

2nr7lh 


(1.8) 


The  method  for  calculation  of  h,  the  convective  heat  transfer  coefficient  to  the  NaK,  is  given  in 
Appendix  A.  The  total  power  transferred  to  the  NaK,  qNaK,  is  given  by  Kreith  and  Black  (1980) 
as 


QNaK  ~ 


T  -T 

1C  1  NaK 

t*n 


(1.9) 


If  the  surface  area  of  all  the  components  is  taken  to  be  27rr7l  (i.e.  neglecting  the  difference  in 
surface  area  of  the  different  layers)  and  equation  1 .9  is  substituted  in  equation  1 .6,  Kmean  is  given 
by 


_  1  JlnfeAQ  Info  As)  InfcAi)  ,  1  1 

",e°"  r7\  k,  *5  K  foj 

The  total  heat  flux  from  the  emitter,  Qin,  is  given  by 


(1.10) 


Qin  ~  Qrad  +  QeE  +  Qcs  •  (1  •  f  1) 

Equilibrium  at  the  emitter  requires  that  QFUEL  and  Qin  are  equal.  The  heat  flux  into  the  collector, 
Qout  is  given  by 


Qout  ~  Qrad  +  QeC  +  Qcs  •  (1-12) 

Equilibrium  at  the  collector  requires  that  Qout  and  QNaK  are  equal. 

It  can  be  seen  from  equations  1.3,  1.11  and  1.12  that  at  equilibrium  the  difference  between  QTISA 
and  QNaK  is  equal  to  the  difference  between  QeE  and  QeC,  the  output  power  flux,  JoutVout. 


In  the  Ya-21u  data  there  are  reliable  values  for  ^NaK  at  the  inlet  and  outlet  of  the  reactor  core,  as  a 
function  of  qFUEL.  A  plot  of  this  data  and  the  analytical  expression  used  in  the  codes  are  shown  in 
figure  Cl.  This  data  was  taken  from  run  8  of  Ya-21u.  The  codes  take  TnoK  as  the  average  of  the 
inlet  and  outlet  temperatures,  a  suitable  value  for  a  position  halfway  along  the  length  of  the  TFE. 
The  use  of  this  TNaK  value  as  a  boundary  condition  has  enabled  us  to  "terminate"  the  heat  transfer 
calculations  at  the  NaK  channel.  Firstly  calculations  are  made  of  the  temperatures  and  heat  fluxes 
in  the  TFE  under  typical  Ya-21u  operating  conditions,  using  EQ1D.  The  experimental  qFUEL/TNaK 
calibrations  for  the  reactor  inlet  and  outlet  are  used  to  fix  the  average  TNaK  for  a  given  TISA 
power.  The  code  is  used  to  find  the  equilibrium  conditions.  These  calculations  give  a  value  for 
Qout  at  that  qFUEL.  Since  this  value  of  Qout  depends  on  Jout  (1.12),  the  calculated  value  of  Jout  must 
be  similar  to  that  when  the  data  was  recorded.  The  value  of  PCs  used  in  the  calculation  is  adjusted 
to  achieve  this  agreement.  Ya-21u  has  been  run  at  relatively  low  output  currents,  where  the 
electron  heat  transfer  term  is  the  order  of  10%  of  the  radiative  term.  EQ1D  is  run  at  values  of 
qFUEL  in  the  range  1.0  -  3.5  kW  (equivalent  to  a  total  reactor  thermal  power  of  37  -  130  kW).  It  is 
then  possible  to  produce  a  calibration  of  Qout  and  ^NaK-  This  calibration  is  valid  under  a  wide 
range  of  operating  conditions,  since  TNafC  is  a  unique  function  of  the  heat  flux  into  the  NaK. 
However,  a  slightly  different  fit  is  used  depending  on  whether  electrode  voltage  drops  are 
included  in  the  calculation.  This  is  because  when  the  calibration  is  generated,  the  calculated  Jout, 
and  hence  the  value  of  Qout,  is  dependent  on  the  magnitude  of  the  voltage  drop.  Due  to  the 
relatively  small  contribution  of  electron  cooling  to  the  total  heat  flux  in  these  calculations,  there 
is  a  very  small  difference  between  these  two  calibrations.  Strictly  speaking,  the  effect  of  NaK 
flow  rate  should  be  accounted  for  (equation  A4),  but  this  thermal  resistance  is  insignificant  in 
comparison  with  the  other  terms  in  equation  1.10.  This  new  Qout  /  TNaK  calibration  is  used  in  the 
more  general  calculations,  performed  using  NAK1D.  It  is  shown  in  figure  C2,  together  with  the 
analytical  expression  used  in  the  codes. 

1 .2.2.1  Effective  emissivitv  of  inter-electrode  gap 

The  effective  emissivity  of  the  inter-electrode  gap,  eeff,  is  defined  by  equation  1.2.  The  value 
used  in  these  models  is  given  by  Nikolaev  (1992)  as 

se#  =0.142 +  (l.26xl0-5r£)  .  (1.13) 

In  the  heat  transfer  codes  all  the  values  of  seff  used  in  a  calculation  can  be  multiplied  by  a  scaling 
factor  entered  by  the  user. 

1.2.3  Cesium  pressure  -  reservoir  temperature  relationship 

The  relationship  between  PCs  and  the  equivalent  cesium  reservoir  temperature,  TR,  needed  to 
produce  this  pressure,  is  given  by  Hatsopoulos  and  Gyftopoulos  (1973)  as 


PCs  =  2.45x1 08r;°'5  exp 


(1.14) 


In  a  TOPAZ-II  reactor,  the  Cs  supply  is  not  a  simple  reservoir  of  the  liquid  metal,  since  this 
would  require  very  precise  temperature  control  to  maintain  a  constant  PCs.  Despite  this,  the 
concept  of  equivalent  reservoir  temperature  is  still  applicable  in  calculating  electrode  work 
functions. 

1.2.4  I-V  characteristics 

The  I-V  characteristic  code  is  central  to  both  of  the  heat  transfer  models.  It  consists  of  separate 
models  for  the  ignited  (discharge)  and  unignited  (diffusion)  modes  of  TFE  operation.  Both  of 
these  models  were  originally  developed  by  McVey  (McVey  and  Rasor,  Dahlberg  et.  al.  1994, 
Rasor  Associates  1991),  as  TECMDL  and  UNIG  respectively.  The  theory  of  both  of  these 
models  will  not  be  described  further  here.  Both  of  these  codes  represent  an  improvement  in 
accuracy  over  the  original  “phenomenological”  model  developed  by  Rasor  (1982).  Details  of  the 
I-V  characteristic  code  are  given  in  Section  3  of  this  report. 

The  present  ignited  mode  code  is  based  on  TECMDL,  with  some  changes  to  suit  the  needs  of 
these  calculations.  The  main  difference  is  that  both  TECMDL  and  UNIG  input  a  value  for  Jout 
and  return  the  corresponding  Vout.  For  the  present  work  it  was  more  convenient  to  input  Vout  and 
calculate  Jout.  In  the  ignited  mode  this  has  changed  the  order  of  some  iterations  in  the  calculation. 
No  modifications  have  been  made  to  the  UNIG  code,  but  an  iterative  routine  has  been  added  to 
precede  it  which  allows  the  calculation  of  the  current  at  the  required  voltage. 

1 .2.4.1  Work  functions 


The  emitter  work  function  (<|)E)  used  in  these  codes  was  taken  from  Hatsopoulos  and  Gyftopoulos 
(1973).  It  is  shown  in  figure  C3,  together  with  the  analytical  expression  used  in  the  codes.  The 
work  function  depends  only  on  TE/TR  and  is  only  applicable  to  Cs  on  W.  It  was  derived  from 
measurements  of  work  function  as  a  function  of  Cs  coverage.  Since  the  thermionic  emission 
process  is  strongly  dependent  on  emitter  work  function,  it  is  inevitable  that  some  differences 
occur  between  experimental  and  calculated  work  function  value,  even  for  a  “pristine”  TFE.  In  the 
heat  transfer  codes  all  the  values  of  tj>E  used  in  a  calculation  can  be  multiplied  by  a  scaling  factor 
entered  by  the  user. 

The  collector  work  function  used  is  based  on  previous  studies  by  other  workers  comparing 
calculations  with  experiment  (McVey  and  Rasor,  Dahlberg  et.  al.  1994).  It  has  been  found  by 
these  workers  to  give  a  better  fit  to  experimental  data  than  measurements  made  on  individual 
materials  outside  of  the  TFE  collector  environment.  It  is  applicable  to  collectors  made  from 
different  metals.  It  is  a  function  of  TC/TR,  and  is  shown  in  figure  C4,  together  with  the  analytical 
expression  used  in  the  codes. 


1 .2.4.2  Voltage  losses  in  electrodes 


The  option  exists  within  all  the  programs  to  allow  for  voltage  drops  due  to  resistive  losses  in  the 
emitter  and  the  collector.  The  derivation  of  the  following  equations  for  voltage  losses  in  the 
electrodes  is  given  in  Appendix  B.  Nikolaev  et.  al.  (1991)  gives  the  resistivities  of  the  emitter 
(pE)  and  the  collector  (pc),  in  a  TOPAZ-II TFE  at  typical  operating  temperatures  as 

p£  =  (4.3xl0"87£)-2.3xl0“5  (Q  cm)  (1.15) 

pc  =  (2.8  x  l(T8rc)-2.0  x  10“6  (Q  cm)  (1.16) 

For  the  active  zone  of  one  electrode,  the  voltage  drop,  Vact  is  given  by 


P  (T)IX 


(1.17) 


where  p(T)  is  the  resistivity  of  the  electrode  at  the  temperature,  T,  of  the  active  zone  (emitter  or 
collector),  I  is  the  current  flowing  from  the  TFE,  X  is  half  of  the  active  zone  length  (15  cm  in  a 
TOPAZ  II  TFE)  and  Axc  is  the  cross  sectional  area  of  the  electrode  (calculated  using  the 
dimensions  shown  in  figure  1.1). 

For  the  inactive  zone  of  one  electrode,  the  voltage  drop,  Vinact,  is  given  by 


Yto,=xf-(p  „r+ipr2)  (i.is) 

where  pH  is  the  resistivity  at  the  hot  end  of  the  inactive  zone  (i.e.  at  the  emitter  or  collector 
temperature),  Y  is  the  length  of  the  inactive  zone  on  one  side  of  the  TFE  (9.25  cm  in  a  TOPAZ  II 
TFE),  and  p  is  the  resistivity  gradient  along  the  inactive  zone,  measured  from  the  hot  to  the  cold 
end.  In  the  calculations  the  temperature  of  the  cold  end  of  both  the  emitter  and  the  collector  is 
taken  to  be  equal  to  the  NaK  temperature.  The  total  voltage  loss  is  given  by  the  sum  of  the  active 
and  inactive  zone  losses  in  the  emitter  and  the  collector. 
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Appendix  A  -  Calculation  of  convective  heat  transfer  to  NaK  coolant 


For  a  cylindrical  geometry,  the  convective  heat  transfer  coefficient,  h,  between  concentric  annuli 
is  given  by  (Burdi  1964) 


h  =  ^[A  +  BPec]  .  (Al) 

kNaK  is  thermal  conductivity  of  the  NaK  coolant,  a  function  of  TNaK.  A  plot  of  kNaK  as  a  function 
of  temperature,  and  the  analytical  expression  used  in  the  codes  is  shown  in  figure  C8  (Paramonov 
and  El-Genk  1994).  D  is  a  characteristic  distance  given  by  (figure  1.1) 

D  =  2(rs-r7)  .  (A2) 

Pe  is  the  Peclet  number  for  the  NaK  coolant,  defined  as  (Kreith  and  Black  1980) 

Pe-Re.Pr  ,  (A3) 

where  Re  is  the  Reynolds  number  and  Pr  is  the  Prandtl  number  of  the  NaK. 

Pe  can  be  expressed  as  (Kreith  and  Black  1980) 


Pe  = 


mDcp 
k  A 

^NaK^NaK 


(A4) 


where  m  is  the  NaK  mass  flow  rate  through  the  TFE,  cp  is  the  NaK  specific  heat  capacity  and 
ANaK  is  the  cross-sectional  area  of  the  NaK  channel,  given  by 

ANaK  =^(r32-f2)  ■  (A5) 

A  plot  of  cp  as  a  function  of  temperature,  and  the  analytical  expression  used  in  the  codes  are 
shown  in  figure  C9  (Paramonov  and  El-Genk  1994). 

Defining  y  as  the  ratio  r8/r7,  the  coefficients  A,  B,  and  C  are  given  by  (Burdi  1964) 

A  =  4.63  +  0.686y  ,  (A6) 

B  =  0.215  +  4.3  x  10"5y  ,  (A7) 


C  =  0.752  +  0.0 1 657y -8.83x1 0"4y 2  . 


(A8) 


Appendix  B  -  Calculation  of  voltage  losses  in  electrodes 


A  schematic  of  the  current  path  through  a  single  TOPAZ-II  TFE  is  shown  in  figure  Bl.  I  is  the 
total  current  through  the  TFE.  There  are  current  connections  at  both  ends  of  the  emitter  and 
collector.  It  has  been  assumed  that  half  of  the  total  current  flows  through  each  electrode 
connection.  The  current  density  is  taken  to  be  uniform  across  the  active  zone  of  the  TFE,  and 
zero  elsewhere. 
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Figure  Bl :  Current  path  through  a  TOPAZ-II  TFE. 

The  assumed  temperature  of  each  electrode  as  a  function  of  distance  along  the  electrode  is  shown 
in  figure  B2. 
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Figure  B2:  Emitter  and  collector  temperature  profiles  used  in  calculations. 


The  variation  of  electrical  resistivity  with  position  along  the  electrodes  is  shown  in  figure  B3, 
together  with  dimensions  used  in  the  following  derivations.  The  resistivity  at  any  point  on  either 
electrode  is  calculated  using  equations  1.15  and  1.16.  Equation  1.15  is  used  for  the  active  zone  of 
the  emitter  (point  Q).  Equation  1.16  is  used  for  the  whole  length  of  the  collector  and  also  the  cold 
ends  of  the  emitter  (point  R).  Resistivities  at  intermediate  points  on  the  inactive  zone  of  the 
emitter  (points  Q  to  R)  are  calculated  using  an  extrapolation  between  the  end  of  the  active  zone 
and  the  cold  end  of  the  emitter.  In  using  the  collector  equation  for  the  cold  end  of  the  emitter,  the 
fact  that  the  emitter  and  collector  are  made  of  different  materials  is  neglected.  However, 
equations  1.15  and  1.16  are  fits  to  experimental  values  of  resistivity  at  typical  emitter  and 
collector  temperatures,  so  it  is  not  valid  to  use  the  emitter  equation  at  the  NaK  temperature.  The 
aim  here  is  to  approximate  the  resistivity  of  the  inactive  zone  rather  than  accurately  model  its 
resistivity. 
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Figure  B3:  Emitter  and  collector  resistivity  profiles  used  in  calculations. 
Active  zone  losses 


At  a  general  point  Xx  along  the  length  of  the  TFE,  the  current  flowing,  i,  (in  either  electrode)  is 
given  by 


I  x 
2~X 


(Bl) 


where  x  is  the  distance  from  the  centre  of  the  TFE  (point  P)  and  X  is  half  the  length  of  the  active 
zone  (15  cm).  The  current  flowing  along  the  electrode  varies  along  the  active  zone,  as  electrons 
are  moving  across  the  gap  along  the  active  zone.  At  either  end  of  the  active  zone,  x  =  X  and  the 
current  =  1/2,  as  shown  in  figure  B 1 . 


Considering  a  small  length  8x  of  one  electrode  at  X,,  the  resistance  8R  of  that  element  is  given 
by 


8R  = 


p(Ty>* 

^xc 


(B2) 


where  Axc  is  the  cross  sectional  area  of  that  electrode.  The  resistivity  of  the  electrode,  p(T)  is 
constant  along  the  active  zone  of  either  electrode,  since  the  temperature  is  assumed  to  be 
constant.  The  voltage  drop,  8V,  due  to  this  resistance  is  given  by  the  product  of  resistance  and 
current  at  that  point, 


5F=Apg> 

2  X  Axc 


(B3) 


The  total  V  drop  from  P  to  Q,  in  either  electrode  is  given  by 


F 


J  2  Y  Axc 


(B4) 


This  gives  the  drop  in  either  active  zone  as 


V„, 


p(r) 


4  A 


IX 


xc 


Inactive  zone  losses 


(B5) 


In  this  region  the  current  is  constant  at  a  value  of  1/2  but  the  resistivity  changes  along  the  length 
due  to  the  change  in  temperature.  Point  Q  (figure  B3)  is  now  considered  as  the  origin  of  x  values, 
for  simplicity.  The  electrode  resistivity  at  Q  is  calculated  from  equation  1.15  or  1.16  as  pQ  and 
that  at  R  is  calculated  from  equation  1 . 1 6  as  pR.  The  value  of  resistivity  at  any  point  X2  is  given 
by 


p(x)=pe  +  px  ,  (B6) 

where  p  is  the  resistivity  gradient  from  Q  to  R  on  the  emitter,  given  by 

P  =  (<0).  (B7) 

Y  is  the  distance  from  the  end  of  the  active  zone  to  the  end  of  the  electrode  (9.25  cm). 

The  resistance  8R  of  element  Sx  on  the  inactive  zone  is  given  by 


The  corresponding  voltage  drop  §V  is  given  by 


8V  = 


1  p(x)8^ 

2  Arc 


(B9) 


By  integrating  equation  B9  and  substituting  equation  B6,  the  total  voltage  drop  Vinact  from  Q  to 
R,  along  the  inactive  zone  of  either  electrode  is  found  to  be 
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Figure  Cl:  NaK  temperature  as  a  fuction  of  TISA  power, 
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Figure  C2:  NaK  temperature  as  a  function  of  TFE  thermal 
output  power  Qout. 


te/t„ 


Figure  C3:  Emitter  work  function,  <j>E,  as  a  function  of 

te/tr. 


T^k  at  the  reactor  outlet  is  fitted  as 


7 —  A  +  Bx  +  Cx ^  +  Dxl  +  Ex  +  Fx^  +  Gx^ , 
x  =  In  {qpuEi) 

A  =  646.08602 
B=  150.47498 
C  =  12.247623 
D=  12.042109 
E  =  10.590804 
F  =  2.4922993 
G  =  0.1794822. 

TNaK  at  the  reactor  inlet  is  fitted  as 

TnoK  ~  A  +  BqFUEL  +  CxqFUEL  +  DqlFUEL  +  Ex/ qFUEL 

A  =  314.80292 
B=  179.19923 
C  = -197.56921 
D=  110.60518 
E  = -7.0386418  x  10'!6. 


TnoK  -  (A  +  B  In  Om„  )  , 


With  electrode 

losses 

Without  electrode 

losses 

A  2.1088286x  10" 

-3.3650990  x  10“* 

B  2.1084566  x  10" 

-3.3642567  x  10*4 

4  +  Cln(7£/Ffi)+E(ln(r£/^)j; 
l  +  B\n(TE/TR)+D(\n(TE/TR)f 

A  =  1.9480958 
B  = -1.2642181 
C  =  -2.8772769 
D  =  0.48526772 
E=  1.4189757. 


$c  =  A  +  B(Tc/TRY  +  cjTyTR 

A  =  214.99825 
B  =  20.424012 
C  =  -176.92743 
D  = -153.3579. 
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Figure  C4:  Collector  work  function,  <j)c,  as  a  function  of 

tc/tr. 


kA  =  A  +  BT  +  CT2  +DT3  +ET* 

A  =  0.749 
B  =  -1.78  x  10’3 
C=  1.75  x  10'6 
D  =  -7.78  x  10'10 
E  =  1.31  x  10'13 


Figure  C5:  Thermal  conductivity  of  A1203,  k4,  as  a  function 
of  temperature. 
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Figure  C6:  Thermal  conductivity  of  He,  k5,  as  a  function  of 
temperature. 
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Figure  C7:  Thermal  conductivity  of  stainless  steel,  k6,  as  a 
function  of  temperature. 
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Figure  C8:  Thermal  conductivity  of  NaK,  k7,  as  a  function 
of  temperature. 
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Figure  C9:  Specific  heat  capacity  of  NaK,  cp,  as  a  function 
of  temperature. 
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Figure  CIO:  Thermal  conductivity  of  Cs  vapor,  A,mw,  as  a 
function  of  weighted  mean  temperature  Tmw. 
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Section  2  -  Applications 


2.1.  Introduction 


The  1-D  thermal  models  EQ1D  and  NAK1D  have  been  used  to  analyse  data  recorded  from  the 
TOPAZ  II  reactor  Ya-21u.  As  has  been  previously  mentioned,  this  reactor  has  been  subjected  to 
a  rigourous  series  of  tests,  both  in  Russia  and  the  US.  Leaks  in  several  of  the  TFEs  have  resulted 
in  exposure  of  the  TFEs  to  air  on  several  occasions  when  the  reactor  has  been  cold.  As  a  result  of 
these  extreme  conditions  changes  in  the  performance  of  the  reactor  have  occurred  over  time.  In 
the  following  analysis  the  model  has  been  used  to  examine  the  effects  of  changes  to  various 
operating  parameters  of  the  reactor  to  help  determine  the  effects  of  extended  testing. 


2.2.  Data  analysis 

2.2.1  Experimental  data 

Data  has  been  taken  from  4  runs  of  Ya-21u.  These  are  :- 


Run 

Dates 

3 

12/10/93  -  12/23/93 

6 

4/1/94  -  4/22/94 

7 

8/1/94  -  8/8/94 

8 

8/10/94  -  8/17/94 

All  these  runs  are  before  the  shake  and  vibration  test,  which  was  performed  immediately 
following  run  8.  Data  from  after  the  shake  and  vibration  test  has  not  been  included  in  this  report 
as  it  shows  much  larger  changes  in  reactor  output,  which  may  have  causes  other  than  those 
considered  here.  In  run  6  the  thermal  shield  to  the  radiator  of  the  reactor  was  removed.  The 
purpose  of  this  shield  is  to  increase  the  NaK  coolant  temperature  (TNaK)  in  the  reactor  during 
electrical  testing  in  vacuum.  This  maintains  TNaK  at  a  more  realistic  operating  temperature  during 
thermal  vacuum  testing.  The  shield  was  present  during  runs  3,  7  and  8.  Runs  7  and  8  were 
intended  to  be  a  single  run,  but  an  emergency  condition  forced  a  shutdown  of  the  reactor  at  the 
end  of  run  7.  Hence  data  from  run  7  is  only  available  for  a  few  operating  conditions.  Air 
intrusions  into  the  TFE  inter-electrode  gaps  over  the  time-span  of  this  data  are  known  to  have 
occurred  prior  to  run  3  and  between  runs  3  and  6. 

Figure  2.1  shows  reactor  I-Y  characteristics  recorded  at  nominal  conditions  of  thermal  power  = 
85  kW  and  95  kW,  PCs  =  0.6  Torr  and  0.7  Torr.  The  data  sets  shown  with  solid  connecting  lines 
have  been  fitted  in  the  calculations  described  below.  The  following  conclusions  can  be  drawn. 
There  has  been  a  gradual  decrease  in  output  from  run  3  to  run  8.  In  figure  2.1a  (85  kW,  0.6  Torr) 
where  data  from  run  7  is  present,  it  shows  a  good  agreement  with  run  6  and  a  less  good 
agreement  with  run  8.  This  shows  that  a  significant  change  occurred  between  runs  7  and  8.  This 


can  most  likely  be  attributed  to  the  emergency  shutdown  and  corresponding  rapid  cool  down  of 
the  reactor  at  the  end  of  run  7.  The  data  from  run  6  fits  well  with  the  trends  shown  by  the  other 
runs,  suggesting  that  removal  of  the  thermal  shield  did  not  have  a  large  effect  on  the  reactor 
output  power  under  these  operating  conditions.  The  data  from  run  8  shows  a  significantly  larger 
incremental  decrease  in  output  than  is  seen  between  the  other  runs,  and  also  a  steeper  I-V 
characteristic.  This  may  indicate  that  the  reactor  has  changed  to  operating  in  the  unignited 
(diffusion)  mode  over  some  of  the  voltage  range.  In  addition  the  I-V  characteristics  in  figures 
2.1b  and  2. Id  (  95  kW,  0.6  and  0.7  Torr)  appear  to  show  a  change  in  gradient  at  around  21  V 
which  may  correspond  to  a  transition  between  the  unignited  and  ignited  modes.  It  should  also  be 
noted  that  in  figure  2.1a  where  two  data  sets  are  shown  from  run  3,  there  is  a  significant 
difference  between  the  two  data  sets.  The  cause  of  this  is  unknown,  but  it  may  be  due  to  the  air 
intrusion  which  occurred  prior  to  this  run.  The  data  set  recorded  on  12/20/93  has  been  fitted  in 
the  following  calculations,  as  this  is  more  consistent  with  the  data  available  in  the  runs  where 
only  one  optimization  was  performed. 

2.2.2  Data  fitting 

Preliminary  analysis  of  the  data  shown  in  figure  2.1  indicates  that  using  the  nominal  values  for 
all  parameters  in  NAK1D  produces  current  densities  around  40%  higher  than  the  measured 
values.  Inaccuracies  in  the  TECMDL  code,  on  which  the  I-V  characteristic  part  of  the  current 
model  is  based,  are  not  enough  to  explain  this  discrepancy.  Precise  comparisons  of  TECMDL 
with  experimental  data  are  difficult  due  to  the  difficulty  in  accurately  measuring  experimental 
parameters,  such  as  emitter  temperature  and  emitter  and  collector  work  functions,  which  are 
required  in  the  calculation.  However  existing  comparisons  (McVey  and  Rasor,  Dahlberg  et.  al. 
1994)  indicate  errors  generally  not  greater  than  10%  between  experimental  data  and  TECMDL. 
The  temperatures  of  the  electrodes,  which  are  used  in  calculating  the  I-V  characteristic,  are 
determined  by  the  heat  transfer  part  of  the  model.  The  heat  transfer  codes  use  published  values 
for  the  thermal  conductivities  of  the  necessary  components.  Whilst  there  may  be  some  deviations 
from  these  values,  they  would  be  unlikely  to  produce  the  differences  between  experiment  and 
calculation  observed  for  the  Ya-21u  data.  The  only  “heat  transfer”  term  which  has  been  allowed 
to  vary  in  this  work  is  the  effective  emissivity,  seff,  of  the  inter-electrode  gap.  This  is  a 
combination  of  the  emissivities  of  the  emitter  and  collector,  and  determines  the  net  radiative  heat 
transfer  between  the  emitter  and  the  collector.  The  emissivity  of  a  surface  is  strongly  dependent 
on  the  surface  finish  and  is  likely  to  have  been  affected  by  the  long  test  history  of  Ya-21u.  The 
effective  emissivity  has  been  varied  using  a  multiplying  factor  which  scales  the  default  value. 

In  fitting  data  from  run  6,  the  option  in  NAK1D  to  change  TNaK  was  used.  Experimental  values 
of  TNaK  at  the  reactor  outlet  in  run  6  were  460  K  and  484  K,  at  85  kW  and  95  kW  respectively.  In 
runs  3  and  8  the  TNaK  values  at  these  same  thermal  powers  were  519  K  and  543  K  respectively. 
Hence  the  change  in  TNaK  from  “normal”  conditions  with  the  radiator  shield  in  place  was  -59  K. 

It  has  long  been  suspected  that  the  Cs  throttle  valve  on  Ya-21u  does  not  accurately  supply  the 
nominal  Cs  pressures  indicated  by  the  valve  position.  This  suspicion  is  based  on  the  observed 
reactor  output,  particularly  discrepancies  between  the  theoretical  and  actual  optimum  Cs 


pressures  for  a  given  thermal  power.  Hence  the  Cs  pressure  (PCs)  has  been  varied  in  these 
calculations  to  fit  the  calculated  output  to  the  experimental  data.  The  primary  effect  of  this  PCs 
change  in  terms  of  output  current  density  (Jout)  is  on  the  emitter  work  function  (<j)E),  via  the 
variation  in  Cs  reservoir  temperature  (see  figure  C3).  The  collector  shows  a  much  slower 
variation  in  work  function  (<j>c)  with  PCs  under  typical  operating  conditions  (figure  C4). 

For  a  given  reactor  thermal  power,  the  electrical  output  power  increases  with  PCs  at  low  PCs, 
passes  through  a  maximum  and  then  falls  more  slowly  as  PCs  is  further  increased.  The  value  of 
the  optimum  PCs  increases  as  the  thermal  power  increases.  In  Ya-21u  experimental  values  of  Jout 
at  given  output  volts  (Vout)  have  been  lower  than  expected.  The  value  used  for  PCs  in  the  model 
can  either  be  decreased  or  increased  to  fit  the  experimental  data.  Fitting  the  95  kW,  0.6  Torr  I-V 
characteristic  from  run  3  with  a  PCs  value  higher  than  the  nominal  value  shows  that  a  PCs  value  of 
4.0  Torr  fits  the  value  of  Jout  at  24  V  and  produces  an  I-V  characteristic  of  similar  gradient  to  the 
experimental  data.  However,  if  the  NaK  temperature  used  in  the  model  is  lowered  by  59  K,  Jout  is 
seen  to  decrease  by  around  24%.  This  is  much  larger  than  the  change  which  was  observed  to 
occur  in  run  6.  The  possibility  that  the  Cs  system  is  producing  much  higher  PCs  values  than 
expected  can  effectively  be  ruled  out  by  this  information.  In  the  following  work,  where  PCs 
values  are  changed  to  fit  the  data,  it  has  been  assumed  that  experimental  PCs  values  are  lower 
than  the  nominal  values. 

The  second  variable  in  this  work  has  been  the  emitter  work  function,  <|)E,  since  this  may  have 
changed  from  its  nominal  temperature  dependence  by  repeated  testing.  Variations  in  (j>E  have 
been  simulated  by  varying  a  multiplying  factor  which  directly  scales  the  whole  of  the  emitter 
Rasor  plot  by  a  constant  factor. 

The  final  variable  has  been  the  value  of  seff  in  the  inter-electrode  gap,  which  may  have  also 
changed  from  its  original  value  after  extended  testing.  Variations  in  seff  have  been  simulated  by 
varying  a  multiplying  factor  which  directly  scales  seff  as  given  in  equation  1.13  by  a  constant 
factor. 

Individual  I-V  characteristics  were  fitted  using  the  value  of  Jout  at  the  nearest  experimental  data 
point  to  Vout  =  24  V.  The  remaining  data  points  on  the  I-V  characteristic  were  calculated  by 
changing  only  the  TFE  output  volts.  PCs,  the  <j)E  multiplier  and  the  seff  multiplier  were  left  the 
same  as  that  used  at  24  V.  Table  2.1  shows  the  results  of  fitting  the  indicated  data  sets  by 
modifying  PCs  to  the  value  shown  in  column  4  from  its  nominal  value  (column  3).  The  fits 
produced  are  shown  in  figures  2.2  and  2.3.  <j)E  was  as  shown  in  figure  C3.  Column  5  shows  the 
calculated  value  of  <j>E  at  24  V  at  the  fitted  PCs  value.  <|)E  varies  by  typically  0.04  eV  between  the 
maximum  and  minimum  voltages  on  the  I-V  characteristics.  This  is  due  to  variations  in  the 
emitter  temperature  (TE)  as  a  result  of  changes  in  output  current  and  therefore  electron  cooling. 
Column  6  shows  the  emitter  work  function  which  would  have  been  calculated  using  the  nominal 
operating  PCs.  Column  7  shows  seff  for  the  inter-electrode  gap,  calculated  from  equation  1.13. 


Run 

TISA 

power 

(kW) 

Nominal 

P cs  (Torr) 

PCs  used 
(Torr) 

MeV) 

Nominal  <J>E 
(eV) 

seff 

3 

85 

0.6 

0.37 

3.129 

2.874 

0.1645 

6 

85 

0.6 

0.37 

3.125 

0.1645 

8 

85 

0.6 

0.33 

3.190 

, 

f 

0.1646 

3 

85 

0.7 

0.37 

3.129 

2.832 

0.1645 

3 

95 

0.6 

0.46 

3.161 

2.978 

0.1649 

6 

95 

0.6 

0.45 

3.170 

0.1649 

8 

95 

0.6 

0.43 

3.203 

• 

0.1650 

3 

95 

0.7 

0.46 

3.162 

2.891 

0.1649 

Table  2.1  :  Cs  pressure  values  required  to  fit  experimental  I-V  characteristics. 

Table  2.2  shows  the  results  of  fitting  the  same  data  sets  using  the  nominal  value  of  PCs  and 
modifying  (j)E  from  its  nominal  value.  This  modification  was  achieved  by  multiplying  the  value 
produced  by  the  emitter  Rasor  plot  in  figure  C3  by  the  factor  shown  in  column  4.  Column  5 
shows  the  calculated  value  of  <f)E  at  24  V.  This  again  varies  by  typically  0.015  eV  between  the 
maximum  and  minimum  voltages  fitted.  The  fits  produced  are  shown  in  figures  2.4  and  2.5.  The 
output  current  was  again  fitted  at  24  V. 


Run 

TISA 

power 

(kW) 

Nominal 

PCs  (Torr) 

multiplier 

<MeV) 

Nominal  (j)E 
(eV) 

£eff 

3 

85 

0.6 

1.044 

3.108 

2.874 

0.1644 

6 

85 

0.6 

1.045 

3.114 

0.1644 

8 

85 

0.6 

1.057 

3.178 

0.1645 

3 

85 

0.7 

1.059 

3.104 

2.832 

0.1644 

3 

95 

0.6 

1.024 

3.151 

2.978 

0.1649 

6 

95 

0.6 

1.026 

3.161 

0.1649 

8 

95 

0.6 

1.031 

3.196 

’ 

0.1650 

3 

95 

0.7 

1.039 

3.149 

2.891 

0.1649 

Table  2.2  :  Required  emitter  work  function  multiplier  to  fit  experimental  I-V  characteristics. 

Table  2.3  shows  the  effect  of  fitting  the  data  from  run  3  at  85,  95  kW,  0.6,  0.7  Torr  by  scaling  the 
value  of  seff  by  the  values  shown  in  column  4.  This  produce  the  seff  values  shown  in  column  5. 
PCs  and  <j)E  were  left  at  their  nominal  values.  However,  <j>E  was  different  to  the  nominal  values 
shown  above,  since  the  increase  in  the  radiative  heat  flux  across  the  inter-electrode  gap  results  in 
lower  Te  values  and  therefore  different  values  of  <J>E.  The  fits  produced  are  shown  in  figure  2.6. 
Since  this  method  of  fitting  I-V  characteristics  did  not  produce  as  good  an  agreement  between 


experiment  and  fit  as  the  previous  methods,  only  the  data  sets  produced  in  run  3  were  fitted  in 
this  way.  The  changes  in  gradient  in  the  fitted  I-V  characteristics  are  caused  by  a  transition  from 
the  unignited  mode  to  the  ignited  mode  of  operation. 


Run 

TISA  power 
(kW) 

Nominal  PCs 
(Torr) 

seff  multiplier 

8eff 

fa  (eV) 

3 

85 

0.6 

1.21 

0.1978 

2.791 

3 

85 

0.7 

1.22 

0.1994 

2.746 

3 

95 

0.6 

1.34 

0.2191 

2.793 

3 

95 

0.7 

1.34 

0.2191 

2.752 

Table  2.3:  Required  effective  emissivity  multiplier  to  fit  experimental  I-V  characteristics. 


2.3  Discussion 


The  following  observations  can  be  made:- 

1.  Tables  2.1  and  2.2  show  that  similar  values  of  <j)E  are  obtained  in  fitting  the  experimental 
data,  whether  the  data  is  fitted  by  changing  PCs  or  by  scaling  <j)E.  This  is  to  be  expected,  since  the 
main  effect  of  a  change  in  PCs  in  I-V  characteristic  calculations  is  to  change  the  electrode  work 
functions. 

2.  For  a  given  PCs,  different  factors  are  required  to  scale  (j>E  at  85  and  95  kW.  This  implies  that 
(|)E  has  degraded  from  its  nominal  Rasor  plot,  and  cannot  be  described  by  uniform  scaling  of  that 
plot.  A  smaller  scaling  factor  is  used  at  95  kW,  where  TE  and  therefore  <|>E  are  higher.  This  means 
that  the  dependence  of  (j)E  on  TE  is  not  as  strong  as  that  in  the  emitter  Rasor  plot.  This  may  be 
expected  for  a  surface  which  has  undergone  long-term  changes  during  the  lifetime  of  the  reactor. 

3.  The  same  PCs  cannot  be  used  to  fit  both  85  and  95  kW  data  at  either  0.6  or  0.7  Torr.  Unless 
the  Cs  throttle  valve  position  cannot  be  accurately  reproduced,  the  same  PCs  values  should  be 
produced  at  different  thermal  powers.  If  the  problem  is  not  caused  by  valve  positioning,  this 
observation  again  suggests  that  <j>E  cannot  be  described  by  the  existing  Rasor  plot.  It  is  unlikely 
that  if  there  was  a  problem  with  reproducibility  in  valve  positioning,  that  the  same  trends  would 
have  been  produced  in  the  different  runs  discussed  here.  Even  if  <j)E  is  no  longer  described  by  the 
nominal  Rasor  plot,  it  is  still  also  possible  that  the  Cs  system  is  not  producing  the  nominal  values 
of  PCs.  These  two  effects  cannot  be  separated  using  the  existing  data. 

4.  There  is  less  change  in  the  I-V  characteristic  than  expected,  when  PCs  is  changed  from  0.6  to 
0.7  Torr  at  either  85  or  95  kW.  Jout  at  24  V  can  be  fitted  with  the  same  value  of  PCs  within  0.01 
Torr,  or  with  a  very  similar  value  of  (t>E.  This  suggests  that  either  the  Cs  valve  is  not  producing  a 
significant  change  in  PCs  or  that  <j)E  has  less  dependence  on  PCs  than  in  the  nominal  Rasor  plot. 


5.  The  calculated  I-V  characteristics  are  too  shallow  in  runs  3  and  8,  too  steep  in  run  6,  when 
fitted  by  changing  PCs  or  the  <j>E  multiplier.  Once  Jout  is  fitted  at  24  V,  the  variation  of  (j>E  as  a 
function  of  the  change  in  TE  along  the  I-V  characteristic  is  assumed  to  obey  the  emitter  Rasor 
plot  shown  in  figure  C3.  Based  on  the  points  made  above,  this  may  not  be  a  valid  assumption, 
and  may  explain  the  discrepancies  in  I-V  characteristic  gradient.  The  magnitude  of  the  difference 
in  Te  along  the  I-V  characteristic  is  similar  in  all  runs,  so  another  explanation  is  needed  to 
explain  the  different  trend  seen  in  run  6.  It  may  be  that  since  the  collector  temperature  is  now 
also  ~  60  K  lower  in  run  6  than  in  the  other  runs,  that  discrepancies  from  the  collector  Rasor  plot 
(figure  C4)  are  also  being  observed. 

6.  When  Jout  was  fitted  using  the  seff  multiplier,  the  I-V  characteristics  produced  were  much  too 
steep.  This  suggests  that  the  changes  in  reactor  output  cannot  be  explained  by  changes  in  seff 
alone.  In  any  case,  it  is  unlikely  that  a  change  in  seff  would  occur  without  changes  in  the  other 
parameters  fitted.  Since  the  fits  produced  by  changing  PCs  or  the  <j>E  multiplier  produce  slightly 
too  low  a  gradient  (point  5),  it  may  be  that  some  changes  in  seff  from  the  nominal  value  are  also 
shown  in  this  data. 

7.  The  emitter  Rasor  plot  used  may  itself  have  some  “systematic”  deviation  from  the  true 
cesiated  plot  used.  For  example,  another  Rasor  plot  for  Cs  on  tungsten  in  given  by  Rasor  (1982). 
This  shows  deviations  of  up  to  0.2  eV  from  the  plot  used  in  this  work,  at  values  of  TE/TR  typical 
of  TOPAZ-II  operating  conditions.  An  error  of  this  kind  would  be  shown  by  the  constant  need 
for  a  <f)E  multiplying  factor.  Even  if  this  kind  of  error  was  present  in  these  calculations,  the  values 
of  (j)E  required  to  fit  the  experimental  data  still  cannot  be  described  by  a  single  Rasor  plot. 


2.4  Summary 

The  experimental  data  presented  above  suggests  that  the  nominal  emitter  Rasor  plot  (figure  C3) 
was  not  adhered  to  by  the  TOPAZ-II  TFEs  during  runs  3  to  8.  In  particular,  the  reactor  output  is 
not  as  sensitive  to  changes  in  emitter  temperature  or  cesium  pressure  as  predicted  by  the  Rasor 
plot  for  Cs  on  tungsten.  It  has  long  been  suspected  that  the  Cs  valve  was  not  producing  the 
correct  cesium  pressures  during  these  runs.  This  observation  may  be  explained  entirely  by  the 
deviation  from  the  Rasor  plot  or  there  may  be  discrepancies  in  the  reactor  output  caused  by  both 
work  function  changes  and  errors  in  cesium  pressure  caused  by  the  cesium  throttle  valve.  There 
may  also  be  some  changes  in  the  effective  emissivity  of  the  inter-electrode  gap  from  the  nominal 
values,  but  the  discrepancies  in  reactor  output  cannot  be  explained  by  this  effect  alone. 

The  absence  of  the  radiator  thermal  shield  in  run  6  did  not  produce  a  large  change  in  output 
power.  This  implies  that  the  Cs  pressure  supplied  by  the  reservoir  is  not  significantly  higher  than 
the  nominal  value. 
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Figure  2.1:  Experimental  I-V  characteristics  from  runs  3,  6,  7  and  8 


Thermal  power  =  85  kW,  fitted  by  modifying  Cs  pressure 
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Thermal  power  =  95  kW,  fitted  by  modifying  Cs  pressure 
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Figure  2.4:  Thermal  power  =  85  kW,  fitted  by  scaling  emitter  work  function 
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Thermal  power  =  95  kW,  fitted  by  scaling  emitter  work  function 
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Figure  2.6:  Run  3,  fitted  by  scaling  effective  emissivity  of  inter-electrode  gap 


Section  3  -  Code  Documentation 


3.1  Introduction 


Three  codes  have  been  written  in  order  to  model  the  TOPAZ -II  TFE  output.  They  are  called 
IV ID,  EQ1D  and  NAK1D,  They  are  all  written  in  Fortran.  The  purpose  of  each  of  these  is 
described  in  the  first  section  of  this  report. 


3.1.1  Symbols 


Symbol 

Name  in  code 

Description 

J's 

JEMIDEAL 

Emitter  zero  field  electron  emission  (Richardson  current  density) 

JCIDEAL 

Collector  zero  field  electron  emission  (Richardson  current  density) 

4>e 

PHIEM 

Emitter  work  function 

HUB! 

PHICOL 

Collector  work  function 

^s 

PHIENEW 

Emitter  work  function  modified  due  to  Schottky  effect  (sat.  mode) 

Te 

EMTEMP 

Emitter  temperature 

Tc 

COLTEMP 

Collector  temperature 

mm 

PRES 

Cesium  pressure 

Tr 

TRES 

Cesium  reservoir  temperature 

TnbK 

TNAK 

NaK  coolant  temperature 

Tm 

TMEAN 

Mean  of  emitter  and  collector  temperatures 

Tmw 

TMW 

Weighted  mean  of  emitter  and  collector  temperatures 

TeE 

TEMPEEM 

Emitter  electron  temperature 

TeC 

TEMPEC 

Collector  electron  temperature 

T 

x  em 

TEEMEAN 

Mean  of  emitter  and  collector  electron  temperatures 

4 _ 

DGAP 

Interelectrode  gap  spacing 

DULAMBDA1 

Ratio  of  gap  spacing  to  electron-neutral  mean  free  path 

m 

DULAMBDA 

Ratio  of  gap  spacing  to  electron  mean  free  path 

kB 

KB 

Boltzmann’s  constant  (=  8.617  x  10‘3  eV  K'1) 

A 

AA 

Factor  preceding  Richardson  equation  (=  120.17266  A  cm'"  K'z) 

mm 

JEM 

Electron  emission  from  virtual  emitter  in  obstructed  mode 

h 

JSEM 

Schottky  modified  emitter  current  density  (in  saturated  mode) 

AV 

DELTAV 

Emitter  barrier  height 

Vout 

VOUTWANT 

VOUT 

Required  output  volts 

Calculated  output  volts 

■*out 

JOUT 

Output  current  density 

mmm. 

JUNIG 

Unignited  mode  output  current  density 

JIXRAT 

Ratio  of  emitter  ion  current  density  to  output  current  density  (<0) 

^E^out 

JERAT 

^S^out 

JSRAT 

Jc^out 

JCRAT 

VE 

VEM 

Emitter  sheath-height  (V) 

Vc 

VCOL 

Collector  sheath  height  (V) 

Vd 

VD 

Arc  drop 

Vrad 

VRAD 

Contribution  to  arc  drop  due  to  photon  losses  from  plasma 

R 

R 

Parameter  involved  in  calculating  JE/J0Ut,  Js/Jout 

^nC 

PHINC 

Equilibrium  plasma  potential 

■ 

JNC 

LTE  current  density 

9fuel 

TPOWER 

Thermal  power  to  1  TFE 

Qfuel 

QFUEL 

Heat  flux  to  emitter  from  fuel  or  TISA 

mm 

QIN 

Heat  flux  from  emitter  into  interelectrode  gap 

H 

QOUT 

Heat  flux  into  collector  from  interelectrode  gap 

QNaK 

QNAK 

Heat  flux  from  collector  to  NaK 

Qrad 

QRAD 

Radiative  heat  flux  from  emitter  to  collector 

Qcs 

QV 

Conductive  heat  flux  from  emitter  to  collector 

QeE 

QEE 

Heat  flux  from  emitter  due  to  electron  cooling 

QeC 

QEC 

Heat  flux  to  collector  due  to  electron  heating 

seff 

EPSTFE 

Effective  emissivity  of  emitter 

■mam 

LAMBDA 

Thermal  conductivity  of  Cs  vapour  at  Tmw 

m 

MDOT 

NaK  mass  flow  rate 

EPSTNEW 

Scaling  factor  for  seff 

EMCOR 

Scaling  factor  for  (|)E 

COLCOR 

Scaling  factor  for  <j>c 

3.2  Code  description 
3.2.1  IV1D 

This  is  a  stand-alone  I-V  characteristic  program.  It  involves  no  heat  transfer  calculations.  It 
calculates  Jout  for  the  required  Vout,  given  PCs,  TE  and  Tc.  It  calculates  Jout  for  both  ignited  and 
unignited  mode  operation,  and  gives  the  final  Jout  as  the  higher  of  the  two  values.  A  schematic  of 
the  structure  of  the  code  is  shown  in  figure  3.1.  The  function  of  each  of  the  subroutines  is 
described  below.  The  reader  should  refer  to  the  TECMDL  and  UNIG  documentation  (Rasor  and 
McVey,  Dahlberg  et.  al.  1994,  Rasor  Associates  1991)  for  further  details  of  the  equations  used  in 
these  codes. 


Figure  3.1  :  Structure  of  IV ID  program 
Main  program 

PCs,  Te  and  Tc  are  entered.  ROOT_FIND  is  called  to  calculate  the  equivalent  Cs  reservoir 
temperature  (TR)  from  PCs.  The  PHIEMSUB  and  PHICOLSUB  subroutines  are  called  to  evaluate 
(j)E  and  (|)c.  The  mean  of  the  emitter  and  collector  temperatures,  Tm  is  calculated  and  used  to 
evaluate  the  ratio  dA,ea.  These  values  are  used  later  in  the  I-V  calculations.  The  required  value  of 
Vout  is  entered.  The  user  has  the  option  to  include  the  effect  of  voltage  drops  in  the  electrodes.  If 
they  are  included,  the  user  is  prompted  for  a  NaK  temperature.  This  is  necessary  as  in  the  heat 
transfer  codes  the  ends  of  the  emitter  and  the  collector  are  assumed  to  be  at  the  same  temperature 
as  the  NaK.  This  value  is  not  used  anywhere  else  in  IV ID.  The  default  NaK  temperature  used 
here  is  25  K  cooler  than  the  collector  temperature. 

ROOT  FIND 


Calculates  the  equivalent  Cs  reservoir  temperature  TR  to  the  Cs  pressure  PCs,  using  a  Newton- 
Raphson  iteration. 

PHIEMSUB 


This  subroutine  calculates  <j>E  from  TE  and  TR.  It  also  uses  the  Richardson  equation  (Hatsopoulos 
and  Gyftopoulos  1973)  to  calculate  the  “ideal”  emitter  current  density  J's ,  assuming  no 
additional  voltage  barriers  to  electron  emission. 

PHTCOLSUB 


Performs  equivalent  calculations  for  the  collector. 


JEMVAL 


If  the  emitter  barrier,  AV,  is  greater  than  zero  (i.e.  in  the  obstructed  region  of  the  ignited  mode), 
J's  is  scaled  by  a  factor  exp(-AV/kBTE),  to  give  the  electron  emission  from  the  “virtual  emitter”, 
JE. 

IV 

The  main  “control”  code  for  calculating  Jout  (ignited  or  unignited  mode)  at  the  required  Vout.  It 
calls  ignited  mode  calculations,  and  compares  these  Jout  values  with  those  from  the  unignited 
mode  calculation.  If  the  ignited  mode  does  not  have  a  solution  at  this  required  Vout,  or  if  the 
unignited  mode  gives  a  higher  Jout,  the  unignited  mode  is  used.  The  following  procedure  is  used. 
AV  is  set  to  zero  and  the  OBSTRUCT  subroutine  is  used  to  calculate  Vout  and  Jout  at  the 
transition  point  (where  the  ignited  mode  switches  from  obstructed  to  saturated  mode  operation). 
It  then  sets  AV  =  -1.0  x  1CT 10  (i.e.  saturated  mode,  very  close  to  transition  point)  and  finds  Vout 
and  Jout  there.  If  the  required  Vout  is  between  these  two  values,  it  uses  a  linear  extrapolation 
between  the  two  points  to  find  the  corresponding  Jout.  This  linear  extrapolation  is  used  because 
points  in  this  AV  range  can  be  difficult  to  find  using  an  iterative  method.  The  unignited  mode 
current  density  (Junig)  is  also  calculated  for  the  required  Vout  and  if  this  value  is  larger  it  will  be 
used  as  the  value  for  Jout.  Control  then  returns  to  the  main  program. 

In  most  cases  the  required  Vout  will  not  fall  in  the  range  described  above,  in  which  case  the 
following  procedure  is  used.  Since  there  are  regions  in  which  the  ignited  mode  I-V  characteristic 
does  not  exist,  a  check  must  be  made  of  whether  it  does  exist  at  the  required  Vout.  The  point  in 
the  obstructed  mode  at  which  the  I-V  characteristic  ceases  to  exist  is  characterized  by  a  region 
where,  as  AV  is  increased,  Vout  reaches  a  maximum  and  begin  to  decrease.  This  effect  occurs 
because  the  model  becomes  non-physical  in  this  region.  The  maximum  voltage  which  can  be 
obtained  in  the  ignited  mode  calculation  is  found  by  increasing  AV  in  progressively  smaller 
increments  and  finding  where  the  corresponding  Vout  begin  to  decrease.  If  the  required  Vout  is 
larger  than  this  value,  the  unignited  mode  is  used  and  control  is  returned  to  the  main  program. 

Having  made  these  preliminary  checks,  a  Newton-Raphson  routine  is  used  to  iterate  AV  until  the 
calculated  value  of  Vout  is  equal  to  the  required  value.  If  the  gradient  of  a  plot  of  Vout  against  AV 
is  found  to  be  negative  during  this  process,  this  corresponds  to  going  above  the  maximum  output 
volts  for  the  ignited  mode  and  the  calculations  is  begun  again  at  a  smaller  AV.  If  more  than  20 
loops  are  used  in  an  iteration,  the  starting  AV  is  made  smaller  and  the  sign  changed,  and  a  new 
iteration  is  started.  Once  a  solution  is  found,  it  is  compared  with  the  unignited  mode  current 
density  and  the  larger  value  taken  as  the  output  current  density.  The  emitter  electron  cooling  heat 
flux,  QeE,  is  also  taken  from  the  ignited  or  unignited  mode  calculation  as  appropriate. 

It  should  be  noted  that  the  quantity  AV  has  been  assigned  both  positive  and  negative  values  to 
correspond  to  the  obstructed  and  saturated  modes  respectively.  In  the  obstructed  mode  it  equals 
the  emitter  barrier  height.  It  does  not  have  physical  meaning  in  the  saturated  mode,  but  has  been 
set  up  to  have  an  artificial  relationship  with  the  quantity  JiX/J0Ut,  which  is  varied  to  vary  Vout  in 


the  saturated  mode  calculation.  Given  the  structure  of  the  present  code,  compared  with 
TECMDL,  this  is  a  more  convenient  approach  and  does  not  affect  the  calculated  output  current 
density. 

OBSTRUCT 


Uses  a  Newton-Raphson  routine  to  solve  the  equation 
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for  a  consistent  value  of  JE/J0Ut.  It  uses  subroutines  PARAMS  if  Vc  >  0  or  PVCLTO  if  Vc  <  0,  to 
calculate  the  necessary  parameters,  some  of  which  are  a  function  of  JE/Jout.  When  convergence 
has  occurred,  it  removes  the  electrode  voltage  losses  from  Vout  and  calculates  QeE. 

PARAMS 

Calculates  values  needed  in  OBSTRUCT,  specifically  TeE,  TeC,  Tem,  Vc,  Vrad,  Vd,  VE,  Vout,  R, 
d/le  and  JE/J0Ut.  It  also  considers  whether  local  thermodynamic  equilibrium  (LTE)  considerations 
affect  either  TeC  alone  or  both  TeC  and  TeE.  Subroutine  TSCAL  is  called  first  to  check  the  value 
of  TeC.  If  the  value  of  TeC  has  to  be  modified,  then  subroutine  LTEAV  is  called  to  change  TeE  if 
this  is  also  required. 

PVCLTO 

If  PARAMS  returns  a  value  of  Vc  <  0,  then  this  subroutine  is  used  instead  of  PARAMS  to 
calculate  the  parameters  required  in  OBSTRUCT.  If  LTE  considerations  affect  TeC  and  TeE 
subroutines  TSCVLTO  and  LTEAVVLTO  are  used  in  the  same  way  as  PARAMS  uses  TSCAL 
and  LTEAV. 


SAT 


If  AV  <  0,  then  this  subroutine  is  called  after  OBSTRUCT.  It  solves  the  equation 
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for  a  constant  Js/Jout.  It  uses  PARSAT  to  calculate  the  necessary  values  for  this  iteration,  some  of 
which  are  a  function  of  Js/Jout.  It  then  solves 


Js  =  J's  exp 
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for  a  constant  Jout.  The  value  of  Jjx/Jout  is  given  by 

A_=AK 


When  convergence  has  occurred,  it  subtracts  the  electrode  voltage  losses  from  Vout  and  calculates 
QeE- 


PARSAT 


Calculates  the  values  needed  for  the  first  iteration  in  SAT,  specifically  TeE,  TeC,  Tem,  Vc,  Vrad, 
Vd,  VE,  Vout,  R,  d/A,e  and  Js/J0Uf  testing  this  I-V  characteristic  over  a  wide  range  of  conditions, 
no  case  has  arisen  in  the  saturated  mode  where  Vc  <  0,  so  this  capability  has  not  been  included  in 
this  model.  If  this  situation  does  arise,  a  comment  will  appear  on  the  screen  and  the  program  will 
pause.  This  case  may  only  arise  under  operating  conditions  far  removed  from  those  relevant  to 
TOPAZ-II.  As  in  PARAMS,  tests  for  LTE  restrictions  are  used.  TSCAL  and  LTEAV  are  used  in 
an  analogous  manner  to  their  use  in  PARAMS. 

TSCAL 

In  order  for  LTE  to  be  satisfied  at  the  collector,  TeC  must  be  greater  than  or  equal  to  the  value 
defined  by  LTE  restrictions.  If  this  is  not  the  case  for  the  value  calculated  in  PARAMS,  it  is 
replaced  by  the  LTE  limited  value,  calculated  using  equations  5,  9,  12  and  13  of  McVey  and 
Rasor.  (Following  references  to  equation  numbers  are  all  from  McVey  and  Rasor). 

LTEAV 


If  TeC  was  constrained  by  LTE,  the  value  of  Tem  may  also  be  restricted.  This  in  turn  changes  both 
TeE  and  TeC.  This  calculation  uses  equations  8  and  9  in  calculating  the  new  value  for  Tem,  then 
calculates  the  new  TeC  using  the  same  procedure  as  in  TSCAL.  It  should  be  noted  that  the 
substitution  TeE  =  2Tem-TeC  has  been  used  in  equation  5  for  the  LTE  restricted  value  of  TeC. 

TSCVLTO 


Performs  the  same  function  as  TSCAL,  but  is  used  by  PVCLTO  for  cases  when  Vc  <  0. 
Calculates  the  LTE  restricted  value  for  TeC  using  the  equation  applicable  for  Vc  <  0. 


LTEAVVLTO 


Performs  the  same  function  as  LTEAV,  but  is  used  by  PVCLTO  for  cases  when  Vc  <  0. 
Calculates  TeC  using  the  equation  applicable  for  Vc  <  0.  The  substitution  TeE  =  2Tem-TeC  has 
again  been  used  in  equation  20  for  the  LTE  restricted  value  of  TeC. 

VCOR 

Calculates  the  internal  voltage  drop  in  the  TFE  as  a  result  of  resistive  losses  in  the  emitter  and 
collector.  If  the  user  has  chosen  not  to  include  voltage  drops,  it  returns  a  value  of  zero. 

UNIGNITED 


Performs  iterations  to  find  Junig  at  the  required  Vout,  using  UNIG  to  calculate  Vout  corresponding 
to  a  given  Junig.  It  first  finds  the  value  of  Vout  at  which  Junig  =  0.  If  the  required  Vout  is  greater  than 
this  value,  it  returns  a  current  density  of  zero.  If  this  is  not  the  case  it  then  uses  the  START 
subroutine  to  find  a  value  of  Junig  producing  a  voltage  less  than  Vout.  A  linear  extrapolation 
between  these  points  is  used  to  find  a  suitable  starting  value  of  current  density  for  the  iterative 
procedure.  At  current  densities  above  a  certain  value,  which  depends  on  the  TFE  operating 
parameters,  UNIG  cannot  find  a  corresponding  Vout.  When  this  is  the  case,  it  returns  a  value  of  - 
99.  If  the  required  (Junig,  Vout)  point  is  very  close  to  this  limit,  the  Newton-Raphson  procedure 
may  “hit”  points  in  this  non-physical  region.  This  is  why  the  initial  extrapolation  to  find  a 
starting  point  close  to  the  required  solution  is  used.  If  an  error  does  occur  in  the  Newton-Raphson 
iteration  a  bisection  procedure  is  used  to  find  the  solution.  UNIG  also  returns  a  value  for  QeE 
which  is  used  in  the  heat  transfer  calculations. 

If  START  cannot  find  a  value  of  Junig  large  enough  to  produce  a  voltage  less  than  Vout,  a  linear 
extrapolation  is  used  from  two  current  values  just  less  than  the  value  at  which  UNIG  ceases  to 
find  a  corresponding  voltage.  An  extrapolated  value  for  QeE  is  also  found. 

START 


Looks  for  a  point  on  the  I-V  characteristic  with  output  volts  less  than  the  required  Vout.  It 
increases  Junig  from  zero  using  progressively  smaller  increments  to  find  a  suitable  point.  If  it  finds 
a  potentially  suitable  point  it  moves  to  a  slightly  value  of  Junig  and  calculates  Junig,  Vout  pairs  at  10 
points  up  to  the  value  initially  chosen  as  a  suitable  starting  point.  It  then  checks  that  there  are  no 
discontinuities  in  the  I-V  characteristic  in  this  region,  by  looking  for  anomalously  large  changes 
in  Vout  between  successive  points.  It  has  been  found  in  tests  that  in  the  region  close  to  where 
UNIG  cannot  return  a  value  for  Vout,  some  discontinuities  can  occur.  If  a  discontinuity  is 
detected,  the  new  starting  point  is  set  to  a  Junig  value  just  below  this  point.  If  the  closest  starting 
point  which  can  be  found  still  produces  a  value  for  Vout  greater  than  the  required  value,  a  logical 
variable  is  set  so  that  the  extrapolation  procedure  is  used  in  UNIGNITED. 


IJNIG 


This  subroutine  is  unchanged  from  McVey’s  version  (Rasor  Associates  1991),  except  that  all 
floating  point  variables  have  been  declared  as  double  precision,  to  be  consistent  with  all  the  other 
subroutines.  It  calculates  Vunig  corresponding  to  a  given  Junig,  for  the  unignited  mode  of  TFE 
operation. 

3.2.2  NAK1D 

NAK1D  calculates  TE,  Tc  and  Jout  for  a  TOPAZ-II  TFE  under  physically  reasonable  operating 
conditions.  Once  the  equilibrium  condition  has  been  found,  it  is  possible  to  change  TNaK  and  find 
the  new  equilibrium  TE,  Tc  and  Jout.  The  structure  of  the  code  is  shown  in  figure  3.2.  The  I-V 
characteristic  calculation  is  the  same  as  that  used  in  IV ID  and  will  not  be  described  further.  In 
the  schematic,  the  subroutine  “IV”  is  taken  to  mean  both  that  subroutine  and  all  other  ignited  and 
unignited  mode  subroutines. 


Figure  3.2  :  Structure  of  NAK1D  program 
Main  Program 

NAK1D  has  been  set  up  to  read  its  input  values  from  a  file,  VALUES.DAT.  The  file  contains 
values  for  PCs,  the  TISA  power  for  one  TFE,  qFUEL,  the  NaK  mass  flow  rate,  tn,  starting  values 
for  Te  and  Tc,  the  required  Vout,  the  option  of  including  or  neglecting  electrode  voltage  drops, 
and  scaling  factors  for  (j)E,  <j>c  and  the  effective  emissivity  of  the  interelectrode  gap.  A  scaling 
factor  of  1.00  leaves  that  quantity  unchanged  from  its  default  value.  The  heat  flux  from  the  TISA 
to  the  emitter,  QFUEL,  is  calculated  from  the  TISA  power.  All  values  in  this  file  must  have  a 
decimal  point,  if  integers  are  present  they  will  be  misinterpreted  by  the  program.  The  value  of  m 
is  for  a  single  TFE,  and  is  therefore  1/37  of  the  value  for  the  NaK  flow  rate  through  the  EM 
pump.  Subroutines  EMIT  and  COLIT  are  called  alternately  to  find  the  equilibrium  values  of  TE 
and  Tc  respectively.  After  each  pass  through  this  pair  of  subroutines  the  user  has  the  option  of 
doing  another  pass  or  stopping  at  that  point.  With  each  pass  through  these  subroutines  the  values 
of  Te,  Tc  and  Jout  get  closer  to  a  consistent  value.  After  this  part  of  the  program  is  completed  the 
user  has  the  option  to  change  TNaK  from  the  value  which  has  already  been  calculated.  If  this 
option  is  chosen,  the  subroutines  COLIT2  and  EMIT  are  used  alternately  until  a  new  equilibrium 
set  of  values  for  TE,  Tc  and  Jout  are  found.  Under  conditions  where  the  TFE  is  close  to  the 


transition  between  the  ignited  and  unignited  modes  of  operation,  the  value  of  TE  or  Tc  may  not 
converge.  If  this  occurs,  the  user  should  change  the  starting  values  of  TE  and/or  Tc.  If  the 
calculation  still  does  not  converge,  the  only  solution  is  to  make  a  small  change  to  one  of  the 
operating  parameters,  for  example  the  thermal  power  or  the  output  volts. 

EMIT 


Calculates  the  equilibrium  value  of  TE.  It  uses  IV  to  return  the  value  for  QeE  and  HEATS  to 
return  the  heat  flux  from  the  emitter  into  the  interelectrode  gap,  Qin.  A  Newton-Raphson  iteration 
is  used  to  find  the  value  of  TE  where  Qin  is  equal  to  QFUEL.  The  collector  is  assumed  to  be  at 
equilibrium  in  this  part  of  the  calculation. 

COLIT 


Calculates  the  equilibrium  value  of  Tc.  It  uses  IV  to  return  Jout  and  calculates  the  heat  flux  to  the 
collector,  Qout  as  the  difference  between  QFuel  and  the  output  power  density.  The  emitter  is 
assumed  to  be  at  equilibrium  in  this  part  of  the  calculation.  TNaK  is  calculated  from  Qout  using 
OUTFIT.  The  heat  flux  to  the  NaK,  QNaK  is  calculated  using  KMEAN.  A  Newton-Raphson 
iteration  is  used  to  find  the  Tc  where  Qout  is  equal  to  QNaK. 

COLIT2 


Calculates  the  equilibrium  value  of  Tc  after  TNaK  has  been  changed  by  the  user.  TNaK  is  now 
fixed.  It  uses  IV  to  return  Jout  and  QeE  and  HEATS  to  calculate  Qout.  QNaK  is  calculated  using 
KMEAN.  A  Newton-Raphson  iteration  is  used  to  find  the  Tc  where  Qout  is  equal  to  QNaK.  The 
emitter  is  again  assumed  to  be  at  equilibrium. 

HEATS 


Calculates  the  values  of  the  heat  fluxes  across  the  interelectrode  gap,  Qrad,  QeE,  QeC  and  QCs.  It 
also  calculates  QNaK  using  KMEAN.  Convective  heat  transfer  to  the  NaK  is  also  included  in  this 
term.  The  effective  emissivity,  seff  is  calculated  using  the  current  value  of  TE.  The  value  of  QeE 
returned  by  the  IV  subroutine  is  used  to  calculate  QeC.  The  weighted  mean  temperature  Tmw  is 
calculated  and  used  in  LAMBCAL  to  find  the  thermal  conductivity  of  the  Cs  vapour,  Lmw. 
Finally  values  of  the  total  heat  flux  from  the  emitter  Qin  and  to  the  collector  Qout  are  found. 

LAMBCAL 


Calculates  Lmw  for  a  given  Tmw  and  PCs. 
OUTFIT 


Calculates  TNaK  from  Qout,  using  the  calibration  produced  using  the  program  EQ1D.  A  slightly 
different  calibration  is  used  depending  on  whether  electrode  voltage  drops  are  included. 


KMEAN 


Calculates  the  value  of  QNaK.  The  convective  heat  transfer  to  the  flowing  NaK  is  also  included  in 
this  term.  The  thermal  conductivities  of  the  plasma-sprayed  alumina,  He  gap  and  stainless  steel 
NaK  channel  wall  are  calculated  from  the  relevant  calibrations  of  thermal  conductivity  versus 
temperature.  Since  the  He  gap  provides  the  greatest  thermal  resistance,  the  plasma-sprayed 
alumina  is  taken  to  be  at  Tc,  the  He  gap  is  taken  to  be  at  the  mean  of  Tc  and  TNaK,  and  the 
stainless  steel  wall  is  taken  to  be  at  TNaK- 


3.2.3  EQ1D 

The  structure  of  this  program  is  very  similar  to  NAK1D.  It  is  used  to  produce  the  calibration 
between  Qout  and  Tnok  used  in  NAK1D.  It  uses  experimental  data  of  the  relationship  between 
qFUEL  and  TNaK  to  fix  TNaK  for  a  calculation  at  a  given  qFuEL-  The  value  of  Qout  produced  by  this 
program  is  paired  with  the  TNaIC  value,  as  described  in  section  1  of  this  report.  The  code  does  not 
include  the  option  to  change  Tn3K>  available  in  NAK1D.  The  only  parts  of  the  program  which  are 
different  from  NAK1D  are  the  main  program  and  subroutine  COLIT2.  Subroutines  POWFIT  and 
POWFIT2  are  only  present  in  EQ1D.  Subroutines  COLIT  and  OUTFIT  from  NAK  ID  are  not 
present.  Only  the  parts  which  are  different  are  described  below.  The  I-V  characteristic  calculation 
is  the  same  as  in  IV  ID.  A  schematic  of  the  structure  of  EQ1D  is  shown  in  figure  3.3. 


Figure  3.3  :  Structure  of  EQ1D  program 
Main  Program 

EQ1D  reads  its  input  values  from  the  same  file,  VALUES.DAT,  as  NAK1D.  The  only  difference 
is  that  once  the  value  for  qFUEL  has  been  read  in,  subroutines  POWFIT  and  POWFIT2  are  called 
to  return  values  for  TNaK  at  the  reactor  outlet  and  inlet  respectively.  The  mean  of  these  two  values 
is  calculated  and  used  as  the  fixed  value  for  TNaIC  through  the  rest  of  the  program. 


C0LIT2 


Calculates  the  equilibrium  value  of  Tc.  It  uses  IV  to  return  Jout  and  calculates  Qout  as  the 
difference  between  QFUEL  and  the  output  power  density.  QNaK  is  calculated  using  KMEAN.  A 
Newton-Raphson  iteration  is  used  to  find  the  Tc  where  Qout  is  equal  to  QNaK.  The  emitter  is 
assumed  to  be  at  equilibrium. 

POWFIT 


Calculates  TNaK  at  the  reactor  outlet  from  qFUEL. 
POWFIT2 


Calculates  TNaK  at  the  reactor  inlet  from  qFUEL. 


